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SYSTEM AND METHOD FOR IMAGE 
SENSING AND PROCESSING 



SPECIFICATION 



BACKGROTHVm OF THF. TNVRNTTON 

5 A number of importaat compression standards for still and video 

images employ the discrete cosine transform (DCT). For example. Fig. 1 illustrates a 
standard JPEG algorithm for compressing a still image. In the illustrated algorithm, 
the image is divided into 8x8 pixel blocks of pixel intensity values (e.g., illustrated 
block 102). For each 8x8 block 102, the two-dimensional (2-D) DCT is computed 

10 (step 104). The DCT coefficients are scaled, quantized, and truncated (i.e., rounded 
off) (step 106) to retain only the information that is most important for accurate 
perception by the human eye. For example, because the eye is relatively insensitive to 
high spatial frequencies, and because the largest DCT coefficients are typically those 
representing the lowest spatial frequencies, many of the high-fi:equeacy DCT 

15 coefficients can be roimded to zero in the quantization step 106. The quantized 
coefficients are then entropy encoded — typically using Huffinan encoding — for 
more compact representation the remaining, non-zero coefficients (step 108). 
The above-described compression scheme can, for example, be applied separately to 
different spectral components of a color image — e.g., the red, green and blue pixels in 

20 an RGB image or the limiiiiance«-chrominance values of the image. Because the DCT 
is a linear operation it can be applied separately to any linear combination of RGB 
pixel values. 

The 2-D, NxN point DCT is defined as follows, 

DCT{A}(kl)^ 

n^ nx a\ r 1 r 1 .,^1 ( TV ^ k ' (2n^l)\ ( 7U ' I ' (2m'\'l)\ 
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where: a(0) = ^,a0c) = J^;c = l,2,3,...N-l , (lb) 

I 

and where A denotes the sampled image, « and m denote the spatial sampling indices, 
and k and / denote the spatial frequency indices. Computation of the 8x8 DCT can 
require on the order of (2x8x8)x(8x8) = 8192 multipUcations, although some weU 
5 known algorithms are capable of reducing the number of multiplications by a factor of 
50 or more. Nonetheless, computation of the DCT typically comprises the bulk of the 
I computations required for image compression. Furthermore, although some 
compression technologies — such as JPEG2000 — use wavelet representations 
rather than the DCT, DCT-based technologies are expected to remain in widespread 
1 0 use for the foreseeable future. 

Moreover, in addition to the JPEG standard, which is used for still 
image compression, there are a number of commonly used video compression 
standards - e.g.. Motion JPEG, MPEG(1,2,4), and H.26X - which require 
computation of the DCT of each frame of the video frame sequence. 

15 Currently, in most commercial appUcations, image compression is 

performed by separate digital signal processing circuits which derive DCT 
coefficients based on digitized image data. However, conventional DCT algorithms 
require a substantial amoxmt of computing power and consume a large amount of 
power, which makes such image processing less attractive for devices in which power 

20 conservation is important. Such devices include, for example, mobile camera phones, 
digital cameras, and wireless image sensors for machine health monitoring and 
surveillance. 



SUMMARY OF THE INVENTTON 
25 It is therefore an object of the present invention to provide an image 

sensing and processing system which reduces the number of computations, 
particularly multiplications, required to derive DCT coefficients from image data. 

It is a further object of the present invention to provide such a system 
which reduces the amount of power consumed by the derivation of DCT coefficients-. 
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These and other objects are accomplished by a system which computes 
DCT coefficieats of an image using the Arithmetic Fourier Transform (AFT). The 
AFT method enables computation of the Fourier transform primarily by performing 
additions. Other than pre-scalingofthe pixel data, no multipUcation is required. M 
5 hardware realizations, the greater computational efficiency of the AFT allows savings 
in circuit complexity, size, and power consumption, and also increases processing 
speed. The image is preferably sampled using non-uniformly spaced sensors, 
although non-uniform sampling can also be achieved by interpolation of signals from 
a set of uniformly spaced sensors. . The AFT algorithm can be implemented in either 
10 digital or analog circuitry. The AFT techniques of the present invention, particularly 
the analog implementations, allow vast economies in circuit complexity and power 
consumption. 

In accordance with one aspect of the present invention, incoming light 
is detected by a sensor array comprising at least first and second sensors having first 

15 and second sensor locations, respectively. The first sensor location is proxunate to a 
location of a first extremum of a basis function of a domain transform, a basis 
function having one or more spatial coordinates defined according to the spatial 
coordinate system of the sensor array. The second sensor location is proximate to a 
location of a second extremum of the same basis fimction or a different basis function. 

20 The system includes at least one filter which receives signals from the first and second 
sensors and generates a filtered signal comprising a weighted sum of at least tiie 
signals &om. the first and second sensors. We include the special case of the foregoing 
m which the signal from a single sensor may comprise a filter output. 

In accordance with additional aspect of the present invention, incoming 

25 Ught is detected by a sensor array comprising a plurality of sensors, including at least 
first and second sensors having first and second sensor locations, respectively. The 
incoming light signal has a first value at the first sensor location and second value at 
tiie second sensor location. The system includes an interpolation circuit which 
receives signals from the first and second sensors, tiiese signals representing the first 

30 and second values, respectively, of tiie incoming Ught signal. The interpolation 
circuit interpolates the signals from tiie first and second sensors to generate an 
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interpolated signal. The interpolated signal represents an approximate value of the 
incoming light signal at a location proximate to a first extremmn of at least one basis 
function of a domain transform, the at least one basis function having at least one 
spatial coordinate defined according to the spatial coordinate system of the sensor 
5 array. 

BRffiF DESCRI PTION OF THE DRAWINGS 
Further objects, features, and advantages of the present invention will 
become apparent firom the following detailed description taken in conjunction with 
10 the accompanying drawings showing illustrative embodiments of the present 
invention, in which: 

Fig. 1 is a block diagram illustrating an exemplary prior art image 
processing procedure; 

Fig. 2 is a diagram illustrating data processed in accordance with the 
1 5 present invention; 

Fig. 3 is a diagram and accompanying graphs illustrating an exemplary 
image sampling space and corresponding domain transform basis functions in 
accordance with the present invention; 

Fig. 4 is a graph illustrating error characteristics of an exemplary 
20 system and method for image sensing and processing in accordance with the present 
invention; 

Fig. 5 is a diagram illustrating an exemplary image sampling space in 
accordance with the present invention; 

Fig. 6 is a graph illustrating error characteristics of an exemplary 
25 system and method for image sensing and processing in accordance with the present 
invention; 

Fig. 7 is a graph illustrating error characteristics of an additional 
exemplary system and method for image sensing and processing in accordance with 
the present invention; 
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Fig. 8 is a graph illustrating error characteristics of yet another 
exemplary system and method for image sensing and processing in accordance with 
the present invention; 

Fig. 9 is a diagram illustrating an exemplary image sampling space in 
5 accordance with the present invention; 

Fig. 10 is a diagram illustrating an exemplary sensor array and filter 
arrangement in accordance with the present invention; 

Fig. 1 1 is a flow diagram illustrating an exemplary image sensing and 
processing procedure in accordance with fhe present invention; 
10 Fig. 12 is a flow diagram illustrating an exemplary signal filtering 

procedure for use in the procedure illustrated in Fig. 1 1 ; 

Fig. 13 is a flow diagram illustrating an exemplary image sensing and 
processing procedure in accordance with the present invention; 

Fig. 14 is a flow diagram illustrating an exemplary signal filtering 
15 procedure for use in the procedure illustrated in Fig. 13; 

Fig. 15 is a diagram illustrating an exemplary sensor array and filtering 
circuit in accordance with the present invention; 

Fig. 16 is a flow diagram illustrating an exemplary image sensing and 
processing procedure in accordance with the present invention; 
20 Fig, 17 is a timing diagram associated with Fig. 10, illustrating an 

exemplary timing sequence produced by the clock generator to generate the filtered 
signal S(3,12). 

Fig. 18 is a diagram illustrating an exemplary sensor array and filter 
arrangement in accordance with the present invention; 
25 Throughout the drawings, unless otherwise stated, the same reference 

numerals and characters are used to denote like features, elements^ components, or 
portions of the illustrated embodiments. 



DETAILED DESCRIPTION OF THE INVENTION 
30 An incoming image signal — such as an incoming light pattern firom a 

scene being imaged — can be sampled by an array of sensors such as a charge 
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coupled device (CCD). In accordance with the present invention, the individual 
sensors in the array can be distributed according to a spatial pattern which is 
particularly well suited for mcreasmg tihe efficiency of AFT algorithms. The 
preferred spatial distribution for a 2-D sensor array can be better understood by first 
5 considering the one-dimensional (1-D) case. For example, to find the 1-D AFT that is 
equivalent to an 8-point, 1-D DCT on a xmit interval (0 to 1) of space or time, 12 non- 
uniformly spaced samples should be used. The preferred sampling locations are (0, 
1/4, 2/7, 1/3, 2/5, 1/2, 4/7, 2/3, 3/4, 4/5, 6/7, 1) — although it is to be noted that, if the 
entire signal being sampled includes multiple unit intervals, the first and the last 

10 samplesof each interval are shared with any adjacent unit intervals. In number 
theory. Sections of the form k/j, where k=0,l . . .N-1 and j=l,2. . .N are commonly 
referred to as "Farey fi-action" of order N. It can thus be seen that the above-described 
sampling locations — which provide the preferred set of samples for calculating an 8- 
point DCT based on the corresponding, 12-point AFT — correspond to an even subset 

15 of Farey fi:^tions of order 8 defined as 2k/j, where k=0,l . . A and j=l,2. . .8. 

The above-described signal samples can be used, in conjunction witii a 
fimction known as the Mobius fimction, to compute the AFT of the signal. The 1-D. 
AFT based on the Mobius function is well known; an exemplary derivation of the 
transform can be found in D. W. Tufts, G. Sadasiv, "Arithmetic Fourier Transform and 

20 Adaptive Delta Modulation: a Symbiosis for High Speed Computation," SPIE Vol. 
880 High Speed Computing (1988). The 1-D Mobius function jii(n) is defined as 
follows: 

fii(l) = 1 (2a) 
^i(ii) - (-1)^ if n=(pi) (pz) (P3).-(ps), where pi are different prime numbers (2b) 
25 p.i(n) = 0 if /7^|« for any prime number jp, (2c) 

where the vertical bar notation m\n means that the integer n is divisible by the integer 
m with no remainder. If n can be expressed as the product of s different prime 
numbers, the value of iiiQi) is (-1)**; otherwise, the value is zero. 
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Within a unit interval, the signal A(t) is assumed to be periodic with 
period one. If the signal A{t) is fiirther assumed to be band-limited to a total of iV 
harmonics, its AFT coef&cients are given by: 

(W ) = • Simhtr^) for* = 1,2,3, , AT , (3) 

5 where each S(n,trej) denotes the output of a filter having the following filtering 

function, based on samples A(t^ - — ) which are distributed at locations 

n 

corresponding to respective Farey fractions of the interval 0 to 1 : 

Sin,t,^) = - ^) fo'^ « = 1 A3, N . (4) 



Each of the filter outputs S{n,tr^ is the sum of the respective samples A(t^j- ) 

10 multiplied by the scale factor 1/n, where tre/ is an arbitrary reference time, tref is 
preferably equal to 1 for a unit interval. Each AFT coefficient is the sum of the filter 
outputs of selected filters, weighted by the Mobius function \iiim). 

To process a 2-D input signal such as an inoage or an image portion 
(e.g., a unit sub-image or block), the AFT is extended to two dimensions usuig a 2-D 
1 5 Mobius function \i2 (ri,m) which is defined as follows: 

|X2 (n,/w)= , (5) 

where n and in are positive integers, and \i\{n) is the 1-D Mobius function defined in 
Eqs. 2a, 2b, and 2c. 

Formulae for the 2-D AFT of a zero-mean 2-D input signal A(p,q) — 
20 where p and q are continuous spatial coordinates in a unit range (i.e., a range firom 0 
to 1) - can be represented with respect to any arbitrary reference point (p^^ qref) by a 
2-D Fourier series as follows: 



CO CO 



(P«£ » <l»f ) = -A^w cos(2«- • kp „f + 6',, )cos(2;r • lq„f + ^, ) , (7) 



7 



wo 2005/017816 PCT/US2003/023160 



where (prefi is an arbitrary reference location, preferably (1, 1) for a unit sub- 
image. 

It is assumed that the signal A(p,q) is band-limited to harmonics in 

« 

both spatial dimensions p and q — i.e., the Fourier series coefficients higher than N 
5 are equal to zero. A filter-bank having filters is used to process the image data, 
each filter having the following filtering function: 

S(n,m,p^,q^)^ SS^CP,^ -"^^4?^^ ) (8) 

nnij^t^ ^ n ^ m 



j k j k 
spatial locations ( ^^^qr^ ) of the samples A{p^ — — , q^^ ) processed 



j k 

respective Farey fractions —and — of tiie dimensions of a xmit image block, as is 



where n = 1,2, • . . • AT and m = 1, 2, N. It can be seen firom Eq. (8) that the 

k J k 
—) of the samples A{p^ — q^^ 

10 by the filters are defibaed — relative to the reference location {prefi qref ) — by 

A 

n m 

discussed in fiirther detail below with respect to Fig. 3. 

By replacing the signal A(p,q) in Eq. (8) by its Fourier series given in 
Eqs. (6) and (7), it can be shown that the output of each filter is equal to the sum of a 
1 5 particular set of Fourier series coefficients of A(p^q) : 

«P CP 

J^l kmi Jsueh that IrxwA Oat 

n\J ml* 

A derivation of Eq. (9) is provided in Appendix A attached hereto. 

Based on the assumption that the signal is band-limited, there are no 



more than 



N 












_n_ 







terms that are non-zero, where \ xj denotes the largest integer 



20 which is less than or equal to x. Given Eq. (9) it is possible to prove the following 
relation for the 2-D Fourier series coefficients (a proof is provided in Appendix B 
attached hereto): 



CO CO 



Ou(j'^>9nif) = Ti^Mi(.nhR)'S(mk,nl,p^,q^) k,l = \,2,....N - (10) 



8 



wo 2005/017816 



PCT/US2003/023160 



Furthermore, because of the close relationship between the DCT and 
the Discrete Fourier Transform (DFT), the above-described outputs from the 2-D AFT 
algorithm can be used to calculate the DCT coefficients of a unit sub-image divided 
into NxN uniformly spaced pixels. First, the image sensor array is divided into unit 
5 area blocks of pixels, each block having by definition, a size of 1x1 . The 

photosensitive elements inside each unit area are placed in locations based on a set of 
Farey firactions of the unit block size, to provide the appropriate samples for the filters 
defined in Eq. (8). In order to calculate the filters' outputs, an ^propriate reference 
location (p,^, is chosen. A convenient reference location is at pr^= 1 and qr^= 1 
10 (at a comer of flie unit area). Eq. (8) then becomes: 

S(n.m)^\^Y.Y.M^-^X--) (11) 

where «=1,2,...JV^ and ni=l,2 JV. 

The output of the 2-D AFT is a set of 2-D Fourier series coefficients. 
In order to derive DCT coefficients fi:om the Fourier Series coefficients, an extended 
15 image block.^p,9) is derived by extending the original image block A{p,q) by its own 
mirror image in both directions, as shown in Fig. 2, as follows: 

A(p.2-q)0^p<l,l^q<2, ^^^^ 
A(2-p.2-qJl^p<2.l^q<2 

If the AFT is to be computed firam extended image block X(p,q), rather 
than firom the original block A(p,qX the appropriate filter values are: 

S(...J=llY.t^(2-^^-^, . (13) 

If the extended image block X(p,q) obeys the Nyquist criterion, the 
resulting AFT coefficients are equal to the DCT coefficients within a scale factor — a 
proof of this result is provided in Appendix C attached hereto. On the other hand, if 
the extended image does not satisfy the Nyquist criterion, the 2-D AFT coefficients 
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are only an approximation of the 2-D DCT coefiBcients. This situation is more likely 
to occur for images rich in high-frequency components. However, it is possible to 
improve the approximation using aliasing correction techniques which are discussed 
in further detail below. 
5 In any case, from Eqs. (12) and (13), the respective outputs S(n,m) of 

the filters can be expressed as fpllows: 




where n and m take the values from 1 to N. [x^ denotes the smallest integer which is 
greater than or equal to x. 

10 ' From Eq. 14 it is apparent that there are certain points in the sample 

space that are repeated. As a result, by calculating the DCT rather than the DFT, the 
number of independent points in the 2-D AFT is decreased by nearly one-half For 
example, to calculate an 8x8 point DCT inside the unit sub-image, a set of 12x12 
photosensitive elements per unit area is used. The elements at the edges of the unit 

IS area are shared between adjacent sub-images, thus reducing the effective number of 
points per block to 1 1x1 1 . An exemplary non-uniform sample space 300 is illustrated 
in Fig. 3. In the illustrated example, non-uniformly distributed sample points 348 are 
used for the 2-D AFT calculation. The corresponding effective DCT sample points 
398 are distributed imiformly. 

20 With the image sampled as illustrated in Fig. 3, and using filters whose 

filtering fimctions are defined according to £q. (14), the 2-D AFT coefficients xjt^i can 
be computed as follows: 

00 00 

^kj = ^^1^2(^9^) ' S(mk^nl) forkj = 1,2,....N (15a) 
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^k,o = 2 M ipi) • S(mk, N) for A: = 1,2,,.. JNT 



(15b) 



= • -^C^' "0 for 1 = l,2,-.iV 



(15c) 



(15d) 



where is the meaa value of the image, Xk.i are the 2-D AFT coefficients of the 
extended block image X, xjtjo are the coefficients obtained by calculating the 1«*D AFT 
of the mean values of the rows along the p-axis, and jcqj are the coefficients obtained 
by calculating the 1-D AFT of the mean values of the columns along the q-axis. 

The corresponding DCT coefficients can be computed as follows: 



10 



DCr{A}(0,0) = 8*£[^] 

DCT{AKk,0) = 4^2*Xk,o 
DCT{AKQ,l) = 4^l2^xoJ 



AF=l,2,...JSr-l 



/=l,2,...iSr-l 



(15e) 



(15f) 



(15g) 



15 



20 



DCT{AKk,l)^4^Xkj 



kr-U2,...N-l and 1=1,2 ... NA (15h) 



The above discussion demonstrates that using the 2-D AFT to compute 
the DCT coefficients of an image portion allows the entire computation to be 
performed primarily with addition operations, and with very few multiplication 
operations, thus making the 2-D AFT procedure extremely efficient. The source of 
this increased efficiency can be further understood with reference to Fig, 3. The 
drawing illustrates an exemplary 2-D sample area 300 of a sensor array corresponding 
,to the area of a conventional 8x8 block of pixels 398 arranged in a conventional 
pattern. However, in accordance with the present invention, the illustrated region 300 
has certain preferred locations 348 for use with the above-described 2-D AFT 
technique. The preferred locations 348 correspond to extrema (i.e., maxima) of basis 
functions of the transform being performed. For example, it is well known that the 
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basis functions of a Fourier transform are sine and cosine functions of various 
different frequencies (in the case of a time-varying signal) or wavelengths (in the case 
of a spatially varying signal such as an image). In the case of a cosine transform such 
as a DCT, the basis functions are cosine functions of various frequencies (for time- 

5 varying signals) or wavelengths (for spatially varying signals), as given by Eq. 1 . In 
the exemplary sample area 300 illustrated in Fig. 3, columns 331, 332, 333, 334, 335, 
336, 337, 338, 339, 340, 341, and 342 correspond to the locations of respective 
maxima 301, 302, 303, 304, 305, 306, 307, 308, 309, 310, 311, and 312 of cosine 
basis functions 320, 321, 322, 323, 324, 325, 326, and 327, where the spatial 

10 coordinate q of these basis functions is defined according to the spatial coordinate 
system of either the sensor array or the illustrated region 300. In particular, in the 
illustrated example, the spatial coordinate q of the aforementioned basis functions 
320, 321, 322, 323, 324, 325, 326, and 327 is equal to the horizontal coordinate of the 
sensor array, referenced to the left edge (column 331) of the illustrated region 300. 

15 Similarly, rows 381, 382, 383, 384, 385, 386, 387, 388, 389, 390, 391, and 392 of the 
preferred sample locations 348 correspond to respective extrema 351, 352, 353, 354, 
35.5, 356, 357, 358, 359, 360, 361, and 362 of cosine basis functions 370, 371, 372, 
373, 374, 375, 376, and 377, these basis functions having a vertical spatial coordinate 
p which, similarly to 9, is defined according to the spatial coordinate system of either 

20 the sensor array or the illustrated region 300 thereof. 

The 2-D AFT calculation uses only selected samples such that, for 
each selected sample, the relevant basis function has a value of +1 at the location of 
the sample. Such a sampling pattern allows the simplifying assumption that, when 
computing the AFT coefficients Xk,u the pre-scaled input sensor signals need only be 

25 multipUed by a factor of +1, 0 or -1 — hence the use of the 2-D Mobius function 
ji2(m,n)inEq. (10). 

Fig. 10 illustrates an exemplary portion 1004 of a sensor array 1034, 
along with a filter arrangement 1022 for detecting an incoming signal (e.g., a ligjht 
pattern being received from a scene being imaged) and processing the signal to derive 

30 the respective filter ou^uts S{n,m) in Eq. (14). The sensor array portion 1004 has 
sensors 1002 located in the preferred locations for the AFT calculation, these 
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locations being defined to have vertical and horizontal distances, relative to comer 
pixel 1028, which are equal to various Farey fractions naultiplied by the size 1032 of 
the array portion 1004. Optionally, the filtering can be performed by an analog circuit 
1022 as is illustrated in Fig. 10 or by a digital filter 1502 as is illustrated in Fig. 15. In 
5 either case, column selection operations are preferably performed by a column 

selector 1036 under control of a microprocessor 1018, and the respective filter outputs 
S(n,m) are stored in a memory device such as RAM 1016. 

Regardless of whether an analog filter 1022 or a digital filter 1502 is 
being used to compute the filter outputs S{n,m% the illustrated arrangement can be 

1 0 operated according to the exemplary procedure illustrated in Fig. 1 1 • In the illustrated 
procedure, an incoming signal — e.g., a light pattern firom a scene — ^is received by the 
sensor array 1004 (step 1 102). The incoming signal is detected by the respective 
sensors 1002 of the array 1004 to generate sensor signals (step 1104), and the signals 
are received by the analog or digital filter arrangement 1022 or 1502 (step 1 106). 

15 Respective weighted simis of respective sets of sensor signals are derived to generate 
respective filtered signals (step 1118). For example, a weighted sum of a set of sensor 
signals (e.g., a weighted sum of the respective pixel values 1028, 1029, 1030, and 
1031 from the intersections of the rows 1024 and 1026 with the columns 1044 and 
1046) is derived by the filter 1022 or 1502 to generate a filtered signal S(2,3) (step 

20 1118). 

In the case of an analog filter arrangement 1022, the weighted sums 
derived in steps 1 108 and 1110 can be produced in accordance with the procedure 
illustrated in Fig. 12. In the iUxistrated filtering procedure 1 108 or 1 1 10, the signals 
from the respective sensors are amplified with the appropriate gains to generate 

25 respective ampUfied signals (step 1208). For example, the signal from the first sensor 
1028 in row 1024 and column 1044 is amplified with a first gain to generate a first 
amplified signal (step 1202), the signal from the second sensor 1029 in the row 1024 
and column 1046 is amplified with a second gain to generate a second amplified 
signal (step 1204), etc. The resulting amplified signals are integrated to generate the 

30 filtered signal (step 1206). 
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The operation of the analog filtering circuit 1022 illustrated in Fig. 10 
can be further understood with reference to the tinaing diagram illustrated in Fig. 17. 
First, the microprocessor 1018 determines which filter is to be calculated — i.e., 
selects values for n and m. Given the value m, the appropriate columns and O'^amp are 

» 

5 selected. Then, given the value of n, appropriate 0\t and are selected. An 
exemplary timing cycle for calculating the filter S(3,12) is as follows: 

1. n=3,m=12, 

2. 0^tot=l, OW==<), where i== 1,2,4,5,6,7,12, 

3. Select Column 0 

10 4. 0^2==1, 0^s2=l, other O'q-H) 

5. Transfer char^ to integrator 1010, <C>t=l, O'sj=0 

6. Select Column 1/6 

7. <dSi=1, <dSi=1, other O'sj=0 

8. Transfer charge to integrator 1010, <I>r=l, O^^=0 
15 9. Select Column 1/3 

10. 0\i=U <E)*si=l, other O'q=0 

ft 

11. Transfer charge to integrator 1010, Ot=l, <l>'q=0 

12. Select Column 1/2 

13. 0\i=l, 0^i=l, other 0*1,1=0 

20 14. Transfer charge to integrator 1 0 1 0, Ot=l , O'q=0 

15. Select Column 2/3 

16. 0^i=l, <I>^i=l, other <I>'q-=0 

17. Transfer charge to integrator 1010, <l>t=l, O'sj==0 

18. Select Column 5/6 
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19. 0^i=l, 0*si=l, Other O'sj=0 

I 

20. Transfer charge to integrator 1010, Ot=l, O^sj=0 

21. Sample the mtegrator's output, 0^=1 

22. 0^^ainp=U ^'ainp=0, whete i = 1,2,3,4,5,6,7, 

5 23. Transfer charge to amplifier 1012, <I>s3=0, Ot3=l 

24. Perform AD conversion using ADC 1014 and store the digital value S(3,12) 
in RAM 1016 

25. Reset the integrator 1010 and amplifier 1012 

Once the respective filter outputs S(m,h) are derived, the 2-D AFT 

10 coefficients are derived (step 1 1 12). To derive the AFT coefficients (step 1112), the 
filter outputs are weighted using appropriate values of the Mobius fimction as is 
described above with respect to Eqs. (15a)-(15d) above (step 1114), and the resulting 
weighted signals are added/summed in accordance with Eqs. (15a)-(15d) (step 1116). 

It is to be noted that, if a digital filter 1502 is used, as is illustrated in 

15 Fig. 15, the respective signals firom the sensors 1002 in the array 1004 are preferably 
amplified by amplifiers 1006, and the resulting amplified signals are then received 
(converted to digital values) and processed by the digital filter 1502. Those skilled in 
the art will be familiar with numerous conmiercially available, individually 
programmable, special-purpose digital filters which can easily be programmed by 

20 ordinarily skilled practitioners to perform the mathematical operations described 
above. Because the resolution of the analog-to-digital converter (ADC) 1014 in a 
typical image sensor system is no greater than 12 bits, a 16-bit digital signal processor 
is suitable for use as the digital filter 1502. 

The 2-D AFT is based on the assmnption that the mean intensity value 

25 (a/k/a the "DC" value) of the fiill sub-image, as well as mean value for each row and 
column separately, is zero. If there is a non-zero DC value for a row, column, or the 
entire sub-image, that value is preferably used to derive correction values for 
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adjusting fhe appropriate filter outputs S(n^m). The proper correction amounts for the 
case when the entire sub-image has a non-zero mean B[A], are as follows: 




N 

k 



N 
I 



^ 



fc=l^-JV-l; A(0,Z)=-£[4|^5^(m) 



l^yi'N-l (16a) 



k 




Jff«l 




(16b) 



10 



In addition, correction amounts should be computed if the input signal 
has non-zero mean values in any of the rows or columns (i.e., if x^^o or xo^ are non- 
zero). In the case of non-zero mean values in rows or colurons, it is sufficient to 
correct only Xk,i, where k=l,2..,.N-l and 1=1,2. ...N-1. The correction formula is as 
follows: 



ij 



0,1 



Hi(m) 



(17) 



Hie correction factors A(k,l) and AiocaiUht) are then added to the 
uncorrected 2-D AFT coeffidents Xkji to derive corrected 2-D AFT coefEicients A;(^,l) 
as follows: 



15 iU^O=JCiy + zl(fc,0 



A:=0, l = i;2.....N-lorAr=l,2...N-l,l = 0 (18a) 



Ac(*,i) = Jcw + AlikJ) + AocaAh^ A,l =1,2. . . ..N-1 



(18b) 



As an illustrative example, the 8x8 DCT case will now be considered 
It is not necessary to determine exactly the respective mean values of the entire unit- 
20 area sub-image and of the local rows and columns. Rather, it is sufficient to use 

estimates for these mean values. For the mean value E^A\ of the entire sub-image A, 
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the closest estimate, in terms of least mean-square error, is provided by tiie filter 
output that averages the largest number of points. In the general, NxN case this is 
S(NJ^. Jn the case of an 8x8 DCT, the best estimate of the mean E[A] of the entire 
sub-image ^ is as follows: 



64 




(19) 



For the 8x8 DCT case, the resulting global DC correction values for 
each 2-D AFT coefficient - based on Eqs. 16a and 16b - are provided in Table 1 : 

Table 1 




10 



The correction values for each 2-D AFT coefficient when there are 
non-zero column means and/or row means are provided in Table 2: 



Table 2 



1 

k 


1 


2 


3 


4 


5 


6 


7 


1 


2x0.1+2x1.0 


2xo^+Xi,o 


2X0.3 


2X0.4 


2Xo.5-Xi,o 


2Xo.6-Xi.O 


2Xo.7-Xi,o 


"2~ 


X0,1+2X2,0 


Xo.2+X2,0 


X03 


XO/l 


X0^-X2.0 


X0^-X2/) 


X0.7-X2.0 


3 


2X3,0 


X3.0 


0 


0 


-X3,0 


-X3.0 


-X3.0 
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4 


2X4.0 


X4.0 


0 


0 


-X4,0 


-X4.0 


-X4.O 


5 




-Xo;i+X5jO 


-xoj 


-xo^ 


-X0.5-X5.0 


-Xo,6-X5,o 


-Xo,7-X5,0 


6 


-X0.1+2X6.0 


-X0.2+X6.0 


-Xo,3 


-X0.4 


-Xo,5-X6,o 


-Xo,6-X6,0 


-Xo,7-X6,0 


7 


-Xo,l+2X7.0 


-Xo;i+X7.0 


-X03 


-Xo,4 


-X0.5-X7.0 


-Xo,6-X7,0 


-Xo,7-X7,0 



Given the 2-D AFT coefficients jcy of the extended sub-image X and 
the corrected AFT coefficients Ac(k,l), the DCT coefficients of the sub-image A can 
be calculated. The relations between the respective 8x8 point DCT coefficients 
5 DCT(/r,l) and the corresponding corrected 2-D AFT coefficients Ac(k,l) are provided 
in Table 3: 



Table 3 



1 

k 


0 


1 


2 


3 


4 


5 


6 


7 


0 


Ac(8,8)*8 


Ac(0.1)*4 


Ac(0^)* 
4V2 


A.(0,3) 
*4V2 


Ac(0.4)* 

4S/2 


Ac(0,5)* 
4V2 


Ac(0,6)* 
4V2 


Ac(0,7)* 
4V2 


1 


A,X1,0)* 
4V2 




Ac(1.2)* 
4 


Ac(U) 
*4 


Ac(l,4)* 
4 


Ac(l,5)* 

4 


Ac(l,6)* 

4 


A.(l,7)*4 


2 


Ac(2,0)* 
4V2 


Ac(2,l)*4 


Ac(2^)* 

4 


Ac(2,3) 
*4 


Ac(2,4)* 
4 


Ac(2.5)* 
4 


Ac(2,6)* 
4 


Ac(2,7)*4 


3 


Ac(3,0)* 
4V2 


Ac(3,l)*4 


A«(3.2)* 
4 


Ac(3,3) 
♦4 


Ac(3,4)* 
4 


Ac(3,5)* 
4 


Ae(3,6)* 
4 


A«(3.7)*4 


4 


Ac(4,0)* 
4V2 


Ac(4,l)*4 


A«(4,2)* 
4 


A<<4,3) 
♦4 


A.(4.4)* 
4 


Ac(4,5)* 
4 


Ac(4,6)* 
4 


Ac(4,7)*4 


5 


4V2 


Ac(5,l)*4 


Ae(5.2)* 
4 


Ac(5,3) 
*4 


Ac(5,4)* 
4 


Ac(5,5)* 
4 


Ac(5,6)» 
4 


A.(5,7)*4 


6 


Ao(6,0)* 

4V2 


Ac(6.1)*4 


Ao(6,2)* 
4 


Ac(6,3) 
♦4 


A«(6,4)* 
4 


Ac(6,5)* 
4 


Ac(6.6)* 
4 


A.(6,7)*4 


7 


A.(7,0)* 
4V2 


Ac(7,l)*4 


Ac(7,2)* 
4 


Ac(7,3) 
*4 


Ac(7,4)* 
4 


A«(7,5)* 
4 


A.(7,6)* 
4 


A,(7,7)*4 



10 If the image signal being sampled has high spatial-firequency 

components that are not integer multiples of the unit spatial frequency, aliasing is 
likely to introduce a certain amount of error into the DCT coefficients computed with 
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tiie AFT algorithm. For example, as is illustrated in Fig. 2, there are sub-image 
boundaries 204 within the extended sub-image 202 derived £rom the origmal sub- 
image 102. At each of these boundaries there is likely to be a discontinuity in the first 
derivative of the pixel intensity. The discontinuities tend to increase as the input 
5 signal frequency approaches half the Nyquist sampling frequency. The 

discontinuities also tend to increase as Hie phase of the input signal approaches 7c/2. If 
substantial discontinuities are present, the extended sub-iihage 202 will have 
significant Fourier components at firequencies greater than half the Nyquist frequency. 
It is well known that if the Nyquist criterion is violated due to undersampling of an 

10 image signal or other signal, the high firequency harmonics — i.e., the components 
violating the Nyquist criterion — "fold back" to appear at firequencies below half the 
Nyquist firequency. An image extension, as such shown in Figure 2, does not lead to 
aliasing effects if the ijoput signal is uniformly sampled at steps of l/8th of the unit 
interval. However, due to the non-uniform placement of samples which have 

15 locations based on Farey fi^ctions as is discussed above, aliasing errors may arise in 
DCT coefficients computed based on the AFT. The mean-square-error between 
uniformly sampled input signal values and an approximation of this signal — where 
the approximation is computed by taking the inverse DCT of the AFT-based DCT 
coefficients — provides an indication of the accuracy of the AFT-based procedure. 

20 The amount of error can be significant when processing image signals which have 
substantial hig^-firequency content. Exemplary results for mean-square error as a 
fimction of firequency are illustrated in Fig. 4, which plots, as a function of firequency, 
the mean-square error of the approximation signal obtained by taking the inverse DCT 
of the exemplary DCT coefficients derived by the above-described AFT technique. 

25 The illustrated results demonstrate that the error is greatest in the high-frequency 
components. 

Error caused by undersampling not only directly affects the accuracy 
of filter outputs S(n,m) before any DC correction is applied, but also affects the 
accuracy of the DC correction itself. 
30 An improved estimate for the mean value of the image may be 

obtained from the output of a filter, S, that averages a set of points taken at a spatial 
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frequency fhat is not expected to be present in the spectrum of the extended image X - 
P. Paparao, A. Ghosh, "An Improved Arithmetic Fourier Transform Algorithm," SPIE 
Vol 1347 Optical Information-Processing Systems and Architectures n {1990). As a 
conclusion of the aforementioned paper, the increase in the ord©r of the filter S, used 

5 to calculate the mean value, may improve the mean value estimate. Thus, the mean- 
square error should decrease when filters of order higher than 8 are used to estimate 
the mean value in the above-described 8x8 DCT case. The density and the number of 
photosensitive elements that are averaged increases, when the higher order filters are 
used, so one should choose a filter with the highest realizable order, as limited by the 

10 fabrication technology. A particular fabrication technology limits the smallest 

distance between photosensitive elements, thus limiting the highest realizable filter 
order. Li order to not significantly increase the number of photosensitive elements, 
the ord^ of the filter should be divisible by at least one lower order. If the order of 
the filter is divisible by the lower order, the Farey firactions of the lower order filter 

15 would match to a subset of Farey firactions associated with a higher order filter, thus 
the number of additional photo-sensitive elements would not increase substantially. A 
typical example is the filter S(12,12), where 12 is divisible by 2,3,4,6. A filter of 
order 12 requires no greater number of photosensitive elements than does a filter of 
order 8- However, for a 12*-order filter, the photosensitive elements are preferably 

.20 more densely packed in certain parts of the sub-image, as is illustrated in Fig. 5. In 

general, in an N*-order filter, sample locations may be placed at positions 2j/N, where 
j=0,l,...N/2. The estimated mean-square^-error, where filters of order 12 are used to 
estimate the global and local mean values, is shown in Fig. 6. 

Optionally, photosensitive elements located at the exact Farey firaction 

25 locations can be used' to obtain the sample values for the high-ord^ filter 

computations used to estimate the global and local DC values. Alternatively, or in 
addition, the sample values can be obtained by interpolation of neighboring samples 
using interpolation procedures discussed in fiirfher detail below. Furthemiore, filters 
of order higher than 12 may be used to estimate the DC values. However, there is a 

30 tradeoff associated with using higher-order filters: such filters may entail an increase 
in the number of photosensitive elements and/or a decrease of the spacing between the 
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elements. Moreover, increasing the order of the filters beyond a value of 12 typically 
does not provide significant additional benefit For example. Fig. 7 illustrates the 
mean-square error of a system which uses filters of order 16 to estimate the global and 
local DC values. A visual comparison of Figs. 6 and 7 reveals that the error is 
5 approximately ttie same in botti cases. It is therefore apparent that filters of order 12 
provide a better tradeoff between the number of sample points (or pixel density) and 
the overall accuracy. 

Aliasing errors in the non-DC-corrected filter outputs can be reduced 
by introducing additional pixels into the sensor array, provided that the fabrication 

10 technology allows for a sufficientiy dense pixel distribution. To correct for such 
aliasing, AFT coefficients of order higher than the equivalent imiform sampling 
firequency (i.e. coefficients of order higher than 8 for the 8x8 DCT case) can be used 
to correct the lower-order coefficients. The higher-order coefficients can be obtained 
directly firom supplemental Farey-firaction-spaced sensors, interpolated firom 

15 neighboring pixels, or can be estin^ted as a firaction of the lower-order coefficients — 
methods which are described in fiirther detail below. 

By introducing additional pixels at the precise Farey firaction locations, 
it is possible to calculate the higher-order AFT coefficients exactly, which then may 
be used to correct the lower-order AFT coefficients. For example, if M is the number 

20 of DCT coefficients, and the highest realizable order of the Farey firaction space is N, 
where N >M — Le. N=9,10,l 1,12,. . . for M=8. First, the global and local DC 
corrections A{hl) and id/ocaK^O ^® estimated using the highest order (N) filters as 
described above, and are added to the uncorrected AFT coefficients jcjfe,i as is indicated 
in Eqs. (18a) and (18b), above. The resulting DC-corrected AFT coefficients Ac(k,l), 

25 where k,l = 0, 1 ,2 . . .N- 1 , are used to deteraiine the aliasing-correction values: 

Aaiias(k,l)= 0 k = 0,l,....2M-N; 1 =0,1,....2M-N (20a) 

Aaiias(k,l) = ' Ac(k,2M-l) k = 0,1,. . .,2M-N; 1 = 2M-N+1,. . ..M-1 (20b) 

Aaiias(k,l) = - Ac(2M.k,l) k = 2M-N+l,....M-l; 1 = 0,1,....2M-N (20c) 
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AaiiasOcl) = - Ac(k^M-l) - Ac(2M-k,l) + Ac(2M-k,2M-l) k =2M-N+1 M-1; 

1=2M-N+1 M-l(20d) 

■ 

The correction fonnulae in Eqs. (20a) - (20d) are valid when M is an 
even number — which is usually the case — and 2M is greater then N. Aliasing- 
5 corrected 2-D AFT coefficients Acc(k,l) can then be calculated by adding the above- 
listed aliasing-correction values to the DC-corrected AFT values Ac(k,l): 

Acc(k,l) = Ac(k,l) + Aaiias(k,l) k= 0,1,.... M-1; 1=0,1,..., M-1 (21) 

Fig. 8 illustrates the estimated mean-square error in an exemplary case 
in which higher Farey firaction samples are used to correct for aliasing. In the 

10 illustrated example, a Farey firaction sample space of order 12 (i.e., N=12) has been 
used to provide the pixel values, filters of order 12 have been used to estimate global 
and local DC values and higher-order AFT coefficients (coefficients of order 8,9,10, 
and 1 1) have been used to correct for aUasing as is discussed above with respect to 
Eqs. (20a) - (20d). In this example, the maximum estimated mean-square-error is at 

1 5 firequency (6.5,6.5) and is equal to 0.0273. 

In Figs. 4 and 6-8, the estimated mean-square errors were derived by 
assuming for each firequency point (ft, f2) that the input hnage X is a 2-D cosine with 
firequency (fi/2, yi). The 2-D AFT based 2-D DCT coefficients were calculated for 
such an input, and then an inverse 2-D DCT was calculated to obtain image Y. The 

20 mean-square error between image Y and X was calculated and assigned to the 
frequency point (fi, f2). 

The preferred number of image samples to be used for the AFT 
computation tends to increase substantially as the order of the Farey fraction space is 
increased. For example a total of 46 photosensitive elements per unit interval should 

25 be used when N=12. It may be impractical or expensive to fabricate image sensors 
with such a high pixel density, in which case the higher AFT coefficients are 
preferably estimated by interpolation of adjacent pixels. The Farey sampling points to 
be used for filters of order M,M+1,. . .N-1 can be interpolated either from the available 
set of samples or fix>m the set of samples being processed by a particular filter. 
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preferably the highest order filter N (the 12th-order filter in the example given above). 
In any case, an exemplary interpolation system is discussed in further detail below. 

Jn an additional method for calculating higher-order 2-D AFT 
coefficients (e.g., coefficients of order 8, 9, 10, and 11), the higher-order coefficients 

5 are calculated as a fraction of the neigjiboring higher-order coefficients. Specifically, 
one or more higher-order coefficients are first calculated using exact Farey sampling 
points, and the other higher-order coefficients can be estimated from these exact 
values as follows. Assuming that the image A is band-limited and has no firequency 
components beyond half the Nyquist frequency, the correlation between respective 

10 neighboring, higher-order Fourier series coefficients is typically quite high. 

Moreover, simulations have shown that even Fourier series coefficients (coefficients 
8,10,12 in our example) tend to be highly correlated with each other, and similarly, 
odd Fourier series coefficients (coefficients 9 and 1 1 in our example) tend to be 
highly correlated with each other. Accordingly, if one even, higher-order Fourier 

IS coefficient is known, the other even, higher-order coefficients can be estimated. 

Similarly, if one odd, higher-order Fourier coefficient is known, the other odd, higher 
order coefficients can be estimated. For example, if a Farey firaction space of order 7 
and filters of order 12 are used, filters of order 9 can be used to estimate the odd 
higher-order coefficients, and filters of order 12 can be used to estimate the even 

20 higher-order coefficients. An exemplary sample space suitable for such estimations 
(for aliasing-correction) is illustrated in Fig. 9, where locations of the photosensitive 
elements are defined as Farey firactions 2j/n; j=0,l,2...n-l and n=l,2,3,4,5,6,7,9, and 
12. 

in addition, a single system can combine the above-described 
25 techniques of: (a) adding sensors at higher-order Farey fraction locations, and (b) 
interpolating the values firom existing sensors to estimate the values of the incoming 
signal at the appropriate higher-order locations. For example, as is illustrated in Fig. 
9, if a desired higher-order pixel location 906 is quite close to a lower-order pixel 
location 904, and there is a sensor at the lower-order location 904, it may be 
30 preferable to compute an estimated value for the higher-order pixel 906 by 
interpolation, rather than by placing a sensor at the higiher-order location 906. 
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However, if a desired higher-order location 910 is farther away from the nearest 
lower-order pixels 908 and 9 12, it may be preferable to add an extra sensor to the 
sensor array in the higher-order location 910. 

Fig. 16 provides an overview of an exemplary procedure for image 
5 sensmg and processing in accordance with the present invention. Pixel values 1602 
are processed to calculate the filters S(n,m) according to Eq.(14) above (step 1604). A 
set of uncorrected AFT coefficients xjtj are computed based upon the filter values 
S(n,m) (step 1606). If the entire image and the respective rows and columns have no 
non-zero DC components, no mean value correction is required (step 1608). The 

10 AFT coefficients Xk,i are therefore power-normalized — as is illustrated above in Eqs. 
(15e)-(15h) — to derive the DCT coefficients 1618 (step 1616). If, however, a mean 
value correction is appropriate (step 1608) the mean value correction amounts are 
computed (step 1610) and used to correct the AFT coefficients Xk,i for deriving 
corrected coefficients AcQc,!) (step 1612). If no aliasing correction is required (step 

15 1614), the procedure continues to step 1616. However, if aliasing correction is 
appropriate (step 1614), the aliasing corrections are computed as discussed above 
(step 1620), and used to further correct the DC-corrected AFT coefficients A^kl) for 
deriving aUas-corrected coefficients Acdkl) (step 1622). The DCT coefficients 1618 
are then calculated based on the alias-corrected AFT coefficients Acdkl) (step 1616). 

20 As is discussed above, interpolation of measurements firom 

neighboring sensors in a sensor array can be usefiil for estimating the value of a pixel 
adjacent to the locations of the sensors. For example, referring to the unit area 300 
illustrated in Fig. 3, if the AFT method of the present invention is to be used with a 
conventional sensor array having sensors located in uniformly spaced positions 398, 

25 interpolation can be used to estimate the values of the image at the Farey-fi^tion- 
based locations 348. If the computation is being performed by a digital signal 
processor such as the digital filter 1502 illustrated in Fig. 15, the computation of the 
value at a particular Farey firaction location 345 can, for example, be performed by 
computing an average of the respective values generated by the sensors located at the 

30 nearest uniformly ^aced locations 394, 395, 396, and 397. 
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Fig. 13 illustrates an exemplary procedure for deriving AFT 
coefficients using interpolated pixel values. In the illustrated procedure, an incoming 
image signal is received by a sensor array (step 1302). The sensor array can, for 
example, be a conventional array having sensors with uniformly distributed spatial 

5 locations. The incoming signal is detected by the sensors of the array to generate a 
plurality of sensor signals (1304). The sensor signals are received by an interpolation 
circuit (step 1306) which interpolates the sensor signals (step 1308) — e.g., by 
averaging the signals — to generate a set of interpolated signals which represent the 
pixel values at locations defined by Farey firactions as is discussed above. The 

10 interpolated signals are received by a filter arrangement such as the analog filter 1022 
illustrated in Fig. 10 or the digital filter 1502 illustrated in Fig. 15 (step 1310). The 
filter 1022 or 1502 derives respective weighted sums of respective sets of interpolated 
signals to generate respective filtered signals (step 1316). For example, a weighted 
sum of a first set of interpolated signals is derived to generate a first filtered signal 

15 (step 1312), and a weighted sum of a second set of interpolated signals is derived to 
generate a second filtered signal (step 1314). In the case of an analog filter 
arrangement 1022, the weighted sums derived in steps 13 12 and 1314 can be 
produced in accordance with the procedure illustrated in Fig. 14. In the illustrated 
filtering procedure 1312 or 1314, the interpolated signals firom a particular rows and 

20 columns are amplified with the appropriate gains to generate respective amplified 
signals (step 1408). For example, a first interpolated signal is amplified with a first 
gain to generate a first amplified signal (step 1402), a second interpolated signal is 
amplified with a second gain to generate a second amplified signal (step 1404), etc« 
The resulting amplified signals are integrated to generate the filtered signal (step 

25 1406). 

Once the respective filter outputs S{m,n) are derived, the 2-D AFT 
coefficients are derived (step 1112). To derive the AFT coefficients (step 1112), the 
filter outputs are weighted using appropriate values of a Mobius fimction as is 
described above with respect to Eqs. (15a)-(15d) above (step 1 1 14), and the resulting 
30 weighted signals are added/summed in accordance with Bqs. (15a)-(lSd) (step 1 1 16). 
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Fxirther improvement of computational efficiency can be achieved by 
using an analog circuit to perform the aforementioned interpolation. Fig. 18 
illustrates an exemplary analog int^olation circuit 1804 for interpolating pixel 
values from sensors 1806 of a sensor array portion 1802 to derive additional pixels 
5 1808, 1810, and 1812 (pixels of the row 1814 and column 1816) for use in an AFT 
^ computation in accordance with the present invention. To interpolate pixels 1808 of 
the row 1814, the pixels 1826 of the rows 1818 and 1820 are used. Similarly, to 
interpolate the pixels 1810 of the column 1816, the pixels 1828 of the columns 1822 
and 1 824 are used. Although the pixels of interest are not necessarily equidistant from 

10 their neighboring pixels, they can be approximated as equidistant, which results in a 
0.5% error. Each interpolated pixel value is therefore approximated as the average 
value of the two neighboring pixel values. A special case is the pixel 1812 at the 
location where row 1814 and column 1816 intersect. This pixel value will be 
interpolated as an average value of four neighboring pixels (pixel values 1830 at the 

15 intersections (1818,1822), (1818,1824), (1820,1822), and (1820,1824)). Assuming 
that the pixels of interest are eqmdistant from their nearest neighbors allows a 
minimum number of sampling capacitors to be used. 

An exemplary timing cycle for calculating the filter S(3,12) using the 
interpolation circuit 1804 is provided below: 

20 1, n=3,m=12, 

2. O^intPl, O^nrO, where 1=1,2,4,5,6,7,12, 

3. Select Column 0 

4. <I>*s2=l, 0^82=1, other a>'sj=0 

5. Transfer charge to integrator 1832, Or=l, O^sj^O 
25 6. Select Column 1/6 

7. O^si^l, a>\i=l, other O^sj=0 

8. Transfer charge to integrator 1832, <I)t==l, O^g==0 

9. Select Column 1/3 
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10. 0^si=l, <I>%i=l, other <D^q=0 

11. Transfer charge to integrator 1832, Ot=l, O'g=0 

12. Select Column 1/2 

13. 0^i=l, 0)^1=1, other <E)*q-=0 

14. Transfer charge to integrator 1832, <X>r=l, O'sj^O 

15. Select Column 2/3 

16. <I>'8i=l, 0^1=1, other <D'g-=0 

17. Transfer charge to integrator 1832, Ot=l, O q-=0 

18. Select Colunm 4/S — inteipolation colunan 

19. 0^s2=l, ^^82=1, other <E>'sj=0 - note that values in column 4/5 are sampled 
with gain 2 instead of gain 4 

20. Transfer charge to integrator 1832, ®r=l, ®^8j=0 

21. Select Coliimn 6/7 - interpolation column 

22. ^>^s2=l, 0^s2=l» other O g=0 - note that values in column 6/7 are sampled 
with Ihe gain 2 instead of gain 4 

23. Transfer charge to integrator 1832, 4)t=l, O"q=0 

24. Sample the integrators output, Os3-l 

25. 0^^amp=l, O*amp=0, whcre i = 1,2,3,4,5,6,7, 

26. Transfer charge to amplifier 1834, Os3=0, Ot3=l 

27. Perform AD conversion using ADC 1836 and store the digital value S(3,12) 
in RAM 1838 

28. Reset the integrator and amplifier 
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Table 4 presents a comparison of the computational efficiencies of 
several different methods for computing a 1-D, 8-point DCT, including the AFT 
method of the present invention. The comparison is expressed in terms of the 
respective numbers of various types of operations used to compute the 1-D DCT: 
5 Table 4 





Conventional 
DCT 


Fast DCT in 
[Oien 1977] 


Fast DCT in 
[Narasiniha 1978] 


AFT 


Additions 


56 


26 


NA 


33 


Multiplications 


42 


16 


18 


7 


Total Operations 
(8pt) 1-D DCT 


1820 


282 


350 


82 


Total Operations 
(8x8 ) 2-D DCT 
Normalized to 
AFT 


493 


12 


18 


1 



It can be seen from Table 4 that, in terms of the total number of 
operations, the AFT method of the present invention is approximately 3.4 times as 
efScient as the most efficient prior art method for computing a 1-D DCT. 

10 Furthermore, because the number of total operations in the 2-D case is approximately 
proportional to the square of the number of computations in the 1-D case, the AFT 
method of the present invention is approximately 12 times as efficient as the most 
efficient prior art method for computing a 2-D DCT. In addition, because the 
multiplications in the AFT computation comprise pre-scaling of the respective pixel 

15 intensities by integer values, these multiplications can be readily implemented using 
analog circuits such as the filter 1022 illustrated in Fig. 10. By effectively eliminating 
most of the digital multiplications, such an analog filter 1022 allows the AFT system 
of the present invention to use 73 times fewer computations than the most efficient 
prior art system. 

20 Although the present invention has be^ described in coimection with 

specific exemplary embodiments, it should be understood that various changes, 
substitutions, and alterations can be made to the disclosed embodiments without 
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departing from the spirit and scope of the invention as set forth in the appended 
claims. 
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Appendix A 

This Appendix provides a proof of the following relation: 



J sudi that k such that 
n\} ni/Q 



The outputs of the filters are as follows: 



n 



m 



(A-1) 



(A-2) 



and the Fourier Series extension of the image A is provided by Eqs. (6) and (7), which 
are reproduced as follows: 



CO oo 



15 



*=1 1=1 



(A-3a) 
(A-3b) 



Thus the filters' output formula (Eq. (A-2)) can be written as follows: 

Rearranging the summation order. Equation (A-4) can be written as in (A-5) 
Having the relation (A-6), the filters' outputs become as in (A-7) 



1 ^ „ f M-n^A 



> = 



0 



n\o 

otherwise 



(A-4) 



(A-5) 



(A-6) 



00 00 



20 'S(«-'^i'^.<W)=I2:«.^(P-^«^)= S S«/*0'-^.V) fAT) 

^ Jsuchthatksuckhat K^'U 
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Appendix B 

This Appendix provides a proof of the following relation: 



00 GO 



^k.\iPr^>qref) = Z2^^2('".w) • S{mk,nl,p^,q^^) for k,l = U2,....N (B-l) 



m=l 11=1 



The Kroneker function is deJQned as follows: 



5(n,/«)=l for n=7w, 
S(«,/w)=0 elsewhere. 

The Mobius function \xi and Kroneker function 5 are related as follows: 



(B-2a) 
(B-2b) 



d\(m/n) 



1, 

n 

0, otherwise 



(B-3) 



10 



The values of m and n are positive integers and sununations are carried out over all 
positive integer values of d that exactly divide.the positive iateger m/n. 
In order to prove the relation in Eq. (B-l), Eqs. (8-3) and (9) can be used to derive the 
following relations: 



00 00 



00 00 



( Y 



00 00 



CO oo 



(B-4) 




w=l v=l 



<^w,yiPr^, <lr^)5{w,k) • 8(v,/) 
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Appendix C 



10 



2D DCT and 2DAFT coefficients equality 

Image X is the extended version of the unit area sub-image A (as shown in Figure 1). 
According to the two-dimensional case of the Nyquist reconstruction formula, the 
continuous image X can be represented by its samples as follows. 





sine 


^ 1 


T 


T 






\ J 



(C-1) 



Without losing generality we will assume that the sampling period T is the same in 
both dimensions and equals 1/8* As a result, there are 16x16 samples. Let us assume 
that the image X is periodic with period 2x2 units. Thus, Equation (C-1) can be 
written as follows: 

/ ^ \ r 



/sin ? 



1 ^ 

q-(l6l+m+-)T 

2 



(C-2) 



15 



It can be shoiyn using the inv^e Fourier transfonn and the dual fonn of the Poisson 
formula that the summation of the sine-functions is equal to the right side of Eq. (C- 
3): 



CO 

^ sine 



r 1 ^ 

p-{l6-k+n+-)T 

2 



1 ^ 



.2-nk( T J\ 



(C-3) 



Based on Eq. (C-3), Eq. (C-2) can be written as follows: 



«=Oin=0 L ^ ^ J AOjb^ 



1 ^ ^-^^n 1 « W'T^'J 

— > e 



(C-4) 
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Because X(p, q) is the extended version of the image A(py q), as expressed in (C-5), 

Eq. (C-4) can be rearranged into Eq. (C-6) 

'Aip,q) 0^p<\,Q^q<\, 

A2-P>q) l<.p<2,0^q<l, 

Aipjl-q) 0^p<l,l^q<2, 

A{2-t,2-q) l^p<2,l^q<2 



X{p,q) = \ 



(C-5) 



10 



15 



7 7 



18 

~ y e CO. 



cos'(l7r)co5'| -73- (15 -2m) | . 



16 



(C-6) 



The product tenn of the cosine functioiis is as follows: 



(15 - 2n) 



V 



7C>A:'(2n + l) 
16 



(C-7) 



Replacing the product terms with (C-7) and rearranging the order of the summations 
Equation (C-6) becomes the following: 



^ n=^0 m=0 



(«+i)r,(w+i)r 



C05 



71-^'(2W + 1) 

16 



COS 



7t-/-(2m + l)^ 



16 



(C-8) 



From Eq, (C-8) it can be seen that the (n^) sxmmiation term does not depend on the 
sign of k and 1. Also, according to the dejSnition of the two-dimensional DCT given in 
(C-10), Eq. (C-8) can be written as follows: 

04 „aO JWeO L ^ ^ J 
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COS 



f7C-/.(2w + l)^ 



16 



cOiy(/'7C'^)+ 



cos[^±^^ r^-/- (2m 4-1) 



• 7C • p)cos{l ^n-q) 



(C-9) 



The dejSnitioii of the two-dimensional DCT is as follows: 

i?C3^(?liO=^a(^).a©J(„+V(m+Vlco^ 

(C-10) 



2A^ 



where: 



10 a(0) = J-l, a(^) = J-|, Ar = l,2,3,...iV-l 



Eq. (C-9) can therefore be written as follows: 



8 M 4V2 

+ Z, C0J(1 • 31 • «f)+ 2^ 2j \ ^ COs{k ' n • p)cOs{l 'K q)(C-l 1) 

In addition, the extended image X{p,q) can be represented by its two- 
dimensional Fourier series: 
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8 8 8 




(C-12) 




where Xkj (fc,/=l,2..8) are 2DAFT coefficieats of the extended image X. The second 
and tiiird terms of Eq. (C-12) are due to the presence of local row and colmnn non- 
zero mean values. The coefficients inside the second terai are calculated as the ID 

5 AFT of tiie row-means; and the coefficients inside the third term are calculated as the 
ID AFT of the column-means. 

Having representations (C-11) and (C-12) of the image X(t,x)j and having 
orthogonal cosine functions in both formulae, we can conclude that the 2DAFT and 
DCT coefficients are equal except for a constant multiplicative factor in each DCT 

10 coefficient: 



DCT{A}(0,0) = S^E[AI 



DCT{AKkfi) = 4yf2^x,^^ 



DCT{A}(0A) = 4^''Xoj 



/ = l,2,3,...iV^-l 



DCT{AKkA) = 4^x,, 



k = l,2,3,...Ar- 1, 1 = l,2,3,...Ar- 1 (C-13) 
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